# ================================================================================== #
# Complaints About Police Misconduct Have Adverse Effects for Black Civilians (PSRM)
# - Script: Appendix A (Background Information)
# - Author: Patrick Kraft (patrickwilli.kraft@uc3m.es)
# ================================================================================== #


# Load raw data, packages, and custom functions ---------------------------

source(here::here("00-func.R"))
load(here("sqf_ccrb.Rdata"))
load(here("out/results.Rdata"))


# Appendix A.1: Contemporaneous Relationship --------------------

m1 <- NULL

## All incidents
m1[[1]] <- lm(ccrb_d ~ sqf_d, data = total)

## Force used
m1[[2]] <- lm(ccrb_d ~ sqf_force_d, data = total)

## Force used + no arrest
m1[[3]] <- lm(ccrb_d ~ sqf_fna_d, data = total)

## add model names
names(m1) <- c("All Incidents","Force Used","Force Used,\n No Arrest")

## create plot
m1_coef <- prePlot(m1) %>% 
  filter(Variables %in% c("sqf_d","sqf_force_d","sqf_fna_d"))
m1_coef$Model <- factor(m1_coef$Model, rev(m1_coef$Model))
p1a <- ggplot(m1_coef, aes(y=Estimate, ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE, x=Model)) +
  plot_default + geom_linerange() + geom_point(size=1) + coord_flip() +
  ggtitle("(A) OLS Coefficients") + ylab("Coefficient") +
  geom_hline(yintercept=0, linetype="dotted") + xlab("Independent Variable (SQF)")

## simulate quantities of interest
m1_pred <- rbind(sim(m1[[1]], iv=data.frame(sqf_d=seq(min(total$sqf_d, na.rm = T), 
                                                      max(total$sqf_d, na.rm=T),
                                                      length.out = 20)), se="newey"), 
                 sim(m1[[2]], iv=data.frame(sqf_force_d=seq(min(total$sqf_force_d, na.rm = T), 
                                                            max(total$sqf_force_d, na.rm=T),
                                                            length.out = 20)), se="newey"), 
                 sim(m1[[3]], iv=data.frame(sqf_fna_d=seq(min(total$sqf_fna_d, na.rm = T),
                                                          max(total$sqf_fna_d, na.rm=T),
                                                          length.out = 20)), se="newey"))
levels(m1_pred$iv) <- c("All Incidents","Force Used","Force Used,\n No Arrest")

## create plot
p1b <- ggplot(m1_pred, aes(x=ivval, y = mean)) +
  plot_default + 
  geom_point(aes(x=sqf_d, y = ccrb_d, shape="All Incidents", col="All Incidents"), 
             size = .5, alpha = .5, data = total) +
  geom_point(aes(x=sqf_force_d, y = ccrb_d, shape="Force Used", col="Force Used"), 
             size = .5, alpha = .5, data = total) +
  geom_point(aes(x=sqf_fna_d, y = ccrb_d, shape="Force Used,\n No Arrest", col="Force Used,\n No Arrest"), 
             size = .5, alpha = .5, data = total) +
  geom_line(aes(lty=iv,col=iv)) + geom_ribbon(aes(ymax=cihi, ymin=cilo, lty=iv), alpha=.2) +
  geom_hline(yintercept = 0, linetype="dotted") +
  geom_vline(xintercept = 0, linetype="dotted") +
  ggtitle("(B) Fitted Regression Lines") +
  xlab("Change in SQF Incidents (x 1000)") +
  ylab("Change in Complaints (x 100)") +
  scale_linetype(name = "SQF Type") +
  scale_color_discrete(name = "SQF Type") +
  scale_shape(name = "SQF Type") +
  theme(legend.position = c(.84,.26))

## save combined plot
ggsave("out/figA1_contemporaneous.png", 
       grid.arrange(p1a,p1b,ncol=2), 
       height=2.5, width=6.5)


# Appendix A.3: Granger Causality Results ---------------------------------

## CCRB complaints
grang1_ccrb <- grangertest(total$ccrb_d ~ total$sqf_d, order=1)
grang2_ccrb <- grangertest(total$ccrb_d ~ total$sqf_d, order=2)
grang3_ccrb <- grangertest(total$ccrb_d ~ total$sqf_d, order=3)
grang4_ccrb <- grangertest(total$ccrb_d ~ total$sqf_d, order=4)
grang5_ccrb <- grangertest(total$ccrb_d ~ total$sqf_d, order=5)

tab_grang_ccrb <- rbind(c(-grang1_ccrb$Df[2], grang1_ccrb$F[2], grang1_ccrb[["Pr(>F)"]][2]), 
                        c(-grang2_ccrb$Df[2], grang2_ccrb$F[2], grang2_ccrb[["Pr(>F)"]][2]), 
                        c(-grang3_ccrb$Df[2], grang3_ccrb$F[2], grang3_ccrb[["Pr(>F)"]][2]), 
                        c(-grang4_ccrb$Df[2], grang4_ccrb$F[2], grang4_ccrb[["Pr(>F)"]][2]), 
                        c(-grang5_ccrb$Df[2], grang5_ccrb$F[2], grang5_ccrb[["Pr(>F)"]][2]))
colnames(tab_grang_ccrb) <- c("Lags","F-Statistic","Pr(>F)")
tab_grang_ccrb <- apply(round(tab_grang_ccrb,3),2,as.character)

print(xtable(tab_grang_ccrb, label = "tab:grang_ccrb", 
             caption = "Granger causality tests of change in SQF incidents predicting subsequent 
             change in CCRB complaints (up to 5 lags)."), 
      table.placement="ht", caption.placement="bottom", 
      include.rownames=F, file="out/tabA1_grang_ccrb.tex")

## SQF incidents
grang1_sqf <- grangertest(total$sqf_d ~ total$ccrb_d, order=1)
grang2_sqf <- grangertest(total$sqf_d ~ total$ccrb_d, order=2)
grang3_sqf <- grangertest(total$sqf_d ~ total$ccrb_d, order=3)
grang4_sqf <- grangertest(total$sqf_d ~ total$ccrb_d, order=4)
grang5_sqf <- grangertest(total$sqf_d ~ total$ccrb_d, order=5)

tab_grang_sqf <- rbind(c(-grang1_sqf$Df[2], grang1_sqf$F[2], grang1_sqf[["Pr(>F)"]][2]), 
                       c(-grang2_sqf$Df[2], grang2_sqf$F[2], grang2_sqf[["Pr(>F)"]][2]), 
                       c(-grang3_sqf$Df[2], grang3_sqf$F[2], grang3_sqf[["Pr(>F)"]][2]), 
                       c(-grang4_sqf$Df[2], grang4_sqf$F[2], grang4_sqf[["Pr(>F)"]][2]), 
                       c(-grang5_sqf$Df[2], grang5_sqf$F[2], grang5_sqf[["Pr(>F)"]][2]))
colnames(tab_grang_sqf) <- c("Lags","F-Statistic","Pr(>F)")
tab_grang_sqf <- apply(round(tab_grang_sqf,3),2,as.character)

print(xtable(tab_grang_sqf, label = "tab:grang_sqf", 
             caption = "Granger causality tests of change in CCRB complaints predicting subsequent 
             change in SQF incidents (up to 5 lags)."), 
      table.placement="ht", caption.placement="bottom", 
      include.rownames=F, file="out/tabA2_grang_sqf.tex")


# Appendix A.4: Equation Balance ------------------------------------------

## Stationarity tests for main endogenous variables
stationarityTab(list("SQF (Total)" = na.omit(total$sqf),
                     "SQF (Differenced)" = na.omit(total$sqf_d),
                     "CCRB (Total)" = na.omit(total$ccrb), 
                     "CCRB (Differenced)" = na.omit(total$ccrb_d),
                     "Arrests (Total)" = na.omit(total$arrests),
                     "Arrests (Differenced)" = na.omit(total$arrests_d))) %>% 
  xtable(caption = "Tests for unit roots vs. stationarity in monthly SQF incidents, CCRB complaints, and arrest time series.", 
         label = "tab:stationarity", align = "rrcccrr") %>% 
  print(file = "out/tabA3_stationarity.tex", include.rownames = F,
        sanitize.colnames.function = function(x){x})

## Stationarity tests for remaining control variables
stationarityTab(list("Media coverage (Total)" = na.omit(controls$media), 
                     "Media coverage (Differenced)" = na.omit(controls$media_d), 
                     "Unemployment Rate (%)" = na.omit(controls$unemp),
                     "Unemployment Rate (Differenced)" = na.omit(controls$unemp_d), 
                     "Overseas Arrivals (x 100,000)" = na.omit(controls$visit_overseas),
                     "Overseas Arrivals (Differenced)" = na.omit(controls$visit_overseas_d), 
                     "Mean Temperature (F)" = na.omit(controls$temp_mean),
                     "Monthly Precipitation (in)" = na.omit(controls$precip_total))) %>% 
  xtable(caption = "Tests for unit roots vs. stationarity in time series used as exogenous control variables.", 
         label = "tab:stationarity_contols", align = "rrcccrr") %>% 
  print(file = "out/tabA4_stationarity_controls.tex", include.rownames = F, 
        sanitize.colnames.function = function(x){x})

## check residuals
m2 %>% 
  map(resid) %>% 
  map(data.frame) %>% 
  map_dfc(slice_tail, n = 139) %>% 
  mutate(Time = row_number()) %>% 
  pivot_longer(-Time, names_to = "variable", values_to = "Residuals") %>% 
  mutate(
    Series = case_when(
      grepl("sqf", variable) ~ "SQF",
      grepl("ccrb", variable) ~ "CCRB",
      grepl("arrests", variable) ~ "Arrests"
    ),
    Series = factor(Series, levels = c("SQF", "CCRB", "Arrests")),
    race = case_when(
      grepl("(10|11|12)", variable) ~ "Latinos (in CCRB & SQF)",
      grepl("(7|8|9)", variable) ~ "Blacks (in CCRB & SQF)",
      grepl("(4|5|6)", variable) ~ "Whites (in CCRB & SQF)",
      grepl("(1|2|3)", variable) ~ "All CCRB Complaints / SQF Incidents"
    ),
    race = factor(race, levels = c("All CCRB Complaints / SQF Incidents", 
                                   "Whites (in CCRB & SQF)", 
                                   "Blacks (in CCRB & SQF)", 
                                   "Latinos (in CCRB & SQF)"))
  ) %>% 
  ggplot(aes(y=Residuals, x=Time, col=Series, lty = Series)) + 
  geom_line() + 
  facet_wrap(~race, ncol = 2, scales = "free_y") +
  plot_default
ggsave("out/figA2_residuals.png",
       height = 3, width = 6.5)

## prepare for plotting
df_resid <- m2 %>% 
  map(resid) %>% 
  map(data.frame) %>% 
  map_dfc(slice_tail, n = 139) %>% 
  mutate(Time = row_number()) %>% 
  pivot_longer(-Time, names_to = "variable", values_to = "Residuals") %>% 
  mutate(
    Series = case_when(
      grepl("sqf", variable) ~ "SQF",
      grepl("ccrb", variable) ~ "CCRB",
      grepl("arrests", variable) ~ "Arrests"
    ),
    race = case_when(
      grepl("(10|11|12)", variable) ~ "Latinos (in CCRB & SQF)",
      grepl("(7|8|9)", variable) ~ "Blacks (in CCRB & SQF)",
      grepl("(4|5|6)", variable) ~ "Whites (in CCRB & SQF)",
      grepl("(1|2|3)", variable) ~ "All CCRB Complaints / SQF Incidents"
    ),
    race = factor(race, levels = c("All CCRB Complaints / SQF Incidents", 
                                   "Whites (in CCRB & SQF)", 
                                   "Blacks (in CCRB & SQF)", 
                                   "Latinos (in CCRB & SQF)"))
  ) %>% 
  filter(Series == "SQF") %>% 
  split(.$variable)

## combine acf plots
acf1 <- ggAcf(df_resid[[1]]$Residuals) +
  ggtitle("All CCRB Complaints / SQF Incidents") +
  plot_default

acf2 <- ggAcf(df_resid[[3]]$Residuals) +
  ggtitle("Whites (in CCRB & SQF)") +
  plot_default

acf3 <- ggAcf(df_resid[[4]]$Residuals) +
  ggtitle("Blacks (in CCRB & SQF)") +
  plot_default

acf4 <- ggAcf(df_resid[[2]]$Residuals) +
  ggtitle("Latinos (in CCRB & SQF)") +
  plot_default

ggsave("out/figA3_acf.png", 
       grid.arrange(acf1, acf2, acf3, acf4), 
       height = 3, width = 5.7)

## combine pacf plots
pacf1 <- ggPacf(df_resid[[1]]$Residuals) +
  ggtitle("All CCRB Complaints / SQF Incidents") +
  plot_default

pacf2 <- ggPacf(df_resid[[3]]$Residuals) +
  ggtitle("Whites (in CCRB & SQF)") +
  plot_default

pacf3 <- ggPacf(df_resid[[4]]$Residuals) +
  ggtitle("Blacks (in CCRB & SQF)") +
  plot_default

pacf4 <- ggPacf(df_resid[[2]]$Residuals) +
  ggtitle("Latinos (in CCRB & SQF)") +
  plot_default

ggsave("out/figA4_pacf.png", 
       grid.arrange(pacf1, pacf2, pacf3, pacf4), 
       height = 3, width = 5.7)

## test for cointegration
gecm1 <- lm(sqf_d ~ lag(sqf) + ccrb_d + lag(ccrb), data = total)
gecm2 <- lm(sqf_d ~ lag(sqf) + ccrb_d + lag(ccrb) + arrests_d, data = total)
gecm3 <- lm(sqf_d ~ lag(sqf) + ccrb_d + lag(ccrb) + arrests_d + 
              media_d + unemp_d + visit_overseas_d + temp_mean + precip_total + jan, 
            data = left_join(total, controls))
tibble(
  `Model` = c("SQF & CCRB", "+ Arrests", "+ Arrests & Controls"),
  `Error Correction` = c(coef(gecm1)["lag(sqf)"],
                         coef(gecm2)["lag(sqf)"],
                         coef(gecm3)["lag(sqf)"]),
  `t Value` = c(summary(gecm1)$coefficients["lag(sqf)","t value"],
                summary(gecm2)$coefficients["lag(sqf)","t value"],
                summary(gecm3)$coefficients["lag(sqf)","t value"]),
  `MacKinnon Value` = c(mkFunc(gecm1$rank-2, length(gecm1$fitted.values)),
                        mkFunc(gecm2$rank-2, length(gecm2$fitted.values)),
                        mkFunc(gecm3$rank-2, length(gecm3$fitted.values))),
  `Cointegration` = rep("No", 3)
) %>% 
  xtable(label = "tab:gecm", align = "lccccc",
         caption = "Test for cointegration based on Generalized Error Correction Models (GECMs).") %>% 
  print(table.placement="ht", caption.placement="bottom", 
        include.rownames=F, file="out/tabA5_gecm.tex")

